function A_pinv = my_inverse_plus(A)
    [m , n] = size(A);
%     number_U = rank(A);
    V = zeros(n , min(m,n)); 
    U = zeros(m , min(m,n));
    V_temp = zeros(n , 1);
    U_temp = zeros(m , 1);
    precision = 1e-10;
    A_fro = norm(A,"fro");
    Error = A_fro^2;
    Error_precision = 1e-12*A_fro;

    for s = 1:n
        U(:,1) = A(:,s);
        u_vector_norm = sqrt(U(:,1)'*U(:,1));
        if u_vector_norm >= precision
            break;
        end
    end
    U(:,1) = U(:,1)/u_vector_norm;
    i = s;
    V(:,1) = A'* U(:,1);
    Error = Error - V(:,1)'*V(:,1);
    k = 2; 
    while Error > Error_precision
        for s = i+1:n
             U_temp = A(:,s); 
             U_temp = MGS_Orth(U(:,1:k - 1), U_temp ,k - 1);
             u_vector_norm = sqrt(U_temp'*U_temp);
             if u_vector_norm >= precision
                 break;
             end
        end
        U_temp = U_temp/u_vector_norm;
        i = s;
        V_temp = A'* U_temp;
        Error = Error - V_temp'*V_temp;
        U(:,k) = U_temp;
        V(:,k) = V_temp;
        k = k + 1;
    end
    U = U(:,1:(k-1));
    V = V(:,1:(k-1));
    V_pinv = my_inverse(V');
    A_pinv = V_pinv*U';
    
end
